GOES 16

Canales en la región del infrarrojo

Canal Longitud de onda Central \(\mu m\) Observaciones Función de peso
7 3.9 Sensible tanto a la radiación del sol como a la emitida por la Tierra Maxímo en superficie
8 6.2 Sensible al vapor de agua Máximo en 425 hPa.
9 6.9 Predomina la absorción por vapor de agua Maxímo cercano a 460 hPa.
10 7.3 Sensible al vapor de agua Máximo en 640 hPa.
11 8.4 Para estudio de microfísica de nubes y plumas volcánicas con ceniza y dióxido de sulfuro Máximo en niveles bajos
12 9.6 Absorción por ozono Máximo en la estratósfera
13 10.3 Ventana radiativa, absorción despresiable en la atmósfera. Cerca del \(\lambda\) de máxima emisión terrestre Máximo en superficie
14 11.2 Ventana radiativa con mayor absorción de vapor de agua
15 12.3 Ventana radiativa con mayor absorción de vapor de agua
16 13.3 Absorción por dióxido de carbono Máximo cerca de superficie

Al parecer las observaciones de GOES vienen en archivos .nc, uno por canal en radiancias.

## Aire claro

## Warning: Removed 8 rows containing missing values (geom_point).

## Warning: Removed 8 rows containing non-finite values (stat_bin).

Con nubes

## Warning: Removed 16 rows containing missing values (geom_point).

## Warning: Removed 16 rows containing non-finite values (stat_bin).

Corrección del BIAS

De acuerdo a la guía del usuario avanzda de GSI, para asimilar radiansas y que estas generen un impacto positivo en el prónosticos, es impresindible corregigir el bias asociado a cada instrumento y canal.

GSI corrige el BIAS sobre la diferencia entre la observación y el background teniendo en cuenta dos posibles fuentes de BIAS:

Ambos set de coeficientes se deben ir actualizando en cada ciclo de asimilación, actualmente esa actualización se hace internamente (siempre que se configuren los argumentos necesarios en el namelist). De acuerdo a la guía los coeficientes para la correcciónn del bias tienen un spin up muy lento y deben ser entranados por semanas o meses. Una alternativa posibles es repetir cada ciclo de asimilación 10-12 veces antes de pasar al siguiente, cada vez actualizando los coeficientes.

Una tercera opción sería probar con los coeficientes que se guardas en GDAS y que usa GFS. No necesariamente es una buena solución ya que son modelos distintos con distinta resolución y dominio. Pero tal vez se puede partir de esos coeficientes para que el spin up sea más corto (en comparación con empezar en 0).

Queda por responder ¿cómo se que los coeficientes ya están bien entrenados y puedo usarlos?

Funciones de peso: http://tropic.ssec.wisc.edu/real-time/amsu/explanation.html

Lectura de coeficientes de GDAS

Maravillosamente los archivos en GDAS son de texto plano y tienen el formanto GSI friendly. Ahora me resulta obvio, pero nunca se sabe!

Entonces, hay un archivo cada 6 horas que se va actualizando.

  • gdas.txxz.abias y gdas.txxz.abias_pc: Según la guía de GSI es “combined satellite angle dependent and mass bias correction coefficient file”. El segundo para canales “pasivos” (que se monitorean).

Al parecer y de acuerdo a las versiones recientes de la guia de usuario, la correción del bias asociado al ángulo y el ¿modelos? ¿backgroud? se hace conjuntamente y por eso los coeficientes se combinan en un solo archivos.

# Es un archivo de texto plano, que emoción!!
satbias_in_file <- "/home/paola.corrales/datosmunin2/nomads.ncdc.noaa.gov/GDAS/201811/20181120/gdas.t00z.abias"

Para que la corrección se haga en conjunto es necesario indicarlo en el namelist: adp_anglebc = .true.

Archivos presentes en la salida de GSI

  • satbias_in: the input coefficients of bias correction for satellite radiance observations.
  • satbias_out: the output coefficients of bias correction for satellite radiance observations after the GSI run.
  • satbias_pc: the input coefficients of bias correction for passive satellite radiance observations
  • satbias_pc.out

Prueba de asimilación

ana_file <- "/home/paola.corrales/datosmunin/EXP/analysis.ensmean"
ana <- ReadNetCDF(ana_file, vars = c(p = "P", "PB", t = "T", qv = "QVAPOR", 
                                 lon = "XLONG", lat = "XLAT")) %>%
    .[, p := p + PB] %>%
    .[, t := tk(t, p)] %>%
    .[, rh := rh(qv, p, t)] %>% 
    .[, td := td(qv, p) + 273.15] %>% 
    .[, ":="(Time = NULL,
             west_east = NULL,
             south_north = NULL,
             qv = NULL,
             PB = NULL)] %>% 
    .[, c("u", "v") := uvmet(ana_file)] %>% 
  .[, file := "ana"]

guess_file <- "/home/paola.corrales/datosmunin/EXP/wrfarw.ensmean"
guess <- ReadNetCDF(guess_file, vars = c(p = "P", "PB", t = "T", qv = "QVAPOR", 
                                 lon = "XLONG", lat = "XLAT")) %>%
    .[, p := p + PB] %>%
    .[, t := tk(t, p)] %>%
    .[, rh := rh(qv, p, t)] %>% 
    .[, td := td(qv, p) + 273.15] %>% 
    .[, ":="(Time = NULL,
             west_east = NULL,
             south_north = NULL,
             qv = NULL,
             PB = NULL)] %>% 
    .[, c("u", "v") := uvmet(ana_file)] %>% 
    .[, file := "guess"]

dt <- rbind(ana, guess)
dt[bottom_top %in% seq(1, 35, 4)] %>% 
  dcast(bottom_top + lon + lat ~ file, value.var = "t") %>% 
  .[, diff_t := ana - guess] %>% 
  ggplot(aes(lon, lat)) +
  geom_point(aes(color = diff_t)) +
  scale_color_divergent() +
  geom_mapa() +
  coord_sf(xlim = c(-78, -50), ylim = c(-42, -19)) +
  facet_wrap(~ bottom_top) +
  theme_minimal()

dt[bottom_top %in% seq(1, 35, 4)] %>% 
  dcast(bottom_top + lon + lat ~ file, value.var = "p") %>% 
  .[, diff_p := ana - guess] %>% 
  ggplot(aes(lon, lat)) +
  geom_point(aes(color = diff_p)) +
  scale_color_divergent() +
  geom_mapa() +
  coord_sf(xlim = c(-78, -50), ylim = c(-42, -19)) +
  facet_wrap(~ bottom_top) +
  theme_minimal()

files <- list.files("analisis/satbias/E7", pattern = "satbias_2", full.names = TRUE)

satbias <- map(files, function(f){
  meta <- unglue::unglue(f, "analisis/satbias/{exp}/satbias_{date}")
  read_satbias(f) %>% 
    as.data.table() %>% 
    .[, ":="(exp = meta[[1]][["exp"]],
             date = ymd_h(meta[[1]][["date"]]))]
  
}) %>% 
  rbindlist()

sensores <- c("amsua_n15",
              "amsua_n18",
              "amsua_metop_a",
              "amsua_aqua",
              "airs_aqua",
              "hirs4_n19",
              "hirs4_metop-a"
              )

satbias[sensor == "amsua_n15" & channel %in% c(1:15)] %>% 
  # melt(measure.vars = paste0("pred", seq_len(1))) %>% 
  ggplot(aes(date, pred1)) +
  geom_line() +
  facet_wrap(~channel, scales = "free_y")

BC_Constant: constant bias correction term BC_Scan_Angle: scan angle bias correction term BC_Cloud_Liquid_Water: CLW bias correction term BC_Lapse_Rate_Squared: square lapse rate bias correction term BC_Lapse_Rate: lapse rate bias correction term BC_Cosine_Latitude_times_Node: node*cos(lat) bias correction term BC_Sine_Latitude: sin(lat) bias correction term BC_Emissivity: emissivity sensitivity bias correction term BC_Fixed_Scan_Position: external scan angle

files <- list.files("analisis/diagfiles/E7", pattern = "asim", full.names = TRUE)
diag <- map(files, function(f){
  meta <- unglue::unglue(f, "analisis/diagfiles/{exp}/asim_{sensor}_{plat}_{date}.ensmean")
  print(f)
  out <- fread(f)
    # .[V10 == 1] %>% 
    
  if (file.size(f) != 0) {
    out[, date := ymd_hms(meta[[1]][["date"]])]
  }
  out
}) %>%
  rbindlist()
## [1] "analisis/diagfiles/E7/asim_airs_aqua_20181120180000.ensmean"
## [1] "analisis/diagfiles/E7/asim_airs_aqua_20181121050000.ensmean"
## [1] "analisis/diagfiles/E7/asim_airs_aqua_20181121170000.ensmean"
## [1] "analisis/diagfiles/E7/asim_airs_aqua_20181121190000.ensmean"
## [1] "analisis/diagfiles/E7/asim_amsua_aqua_20181120180000.ensmean"
## [1] "analisis/diagfiles/E7/asim_amsua_aqua_20181121050000.ensmean"
## [1] "analisis/diagfiles/E7/asim_amsua_aqua_20181121170000.ensmean"
## [1] "analisis/diagfiles/E7/asim_amsua_aqua_20181121190000.ensmean"
## [1] "analisis/diagfiles/E7/asim_amsua_metop-a_20181121010000.ensmean"
## [1] "analisis/diagfiles/E7/asim_amsua_metop-a_20181121130000.ensmean"
## [1] "analisis/diagfiles/E7/asim_amsua_metop-a_20181121140000.ensmean"
## [1] "analisis/diagfiles/E7/asim_amsua_n15_20181120230000.ensmean"
## [1] "analisis/diagfiles/E7/asim_amsua_n15_20181121000000.ensmean"
## [1] "analisis/diagfiles/E7/asim_amsua_n15_20181121100000.ensmean"
## [1] "analisis/diagfiles/E7/asim_amsua_n15_20181121120000.ensmean"
## [1] "analisis/diagfiles/E7/asim_amsua_n15_20181121220000.ensmean"
## [1] "analisis/diagfiles/E7/asim_amsua_n15_20181122000000.ensmean"
## [1] "analisis/diagfiles/E7/asim_amsua_n18_20181121000000.ensmean"
## [1] "analisis/diagfiles/E7/asim_amsua_n18_20181121020000.ensmean"
## [1] "analisis/diagfiles/E7/asim_amsua_n18_20181121120000.ensmean"
## [1] "analisis/diagfiles/E7/asim_amsua_n18_20181121130000.ensmean"
## [1] "analisis/diagfiles/E7/asim_amsua_n18_20181122000000.ensmean"
## [1] "analisis/diagfiles/E7/asim_hirs4_metop-a_20181121010000.ensmean"
## [1] "analisis/diagfiles/E7/asim_hirs4_metop-a_20181121130000.ensmean"
## [1] "analisis/diagfiles/E7/asim_hirs4_metop-a_20181121140000.ensmean"
## [1] "analisis/diagfiles/E7/asim_hirs4_n19_20181120180000.ensmean"
## Warning in fread(f): File 'analisis/diagfiles/E7/
## asim_hirs4_n19_20181120180000.ensmean' has size 0. Returning a NULL data.table.
## [1] "analisis/diagfiles/E7/asim_hirs4_n19_20181120190000.ensmean"
## Warning in fread(f): File 'analisis/diagfiles/E7/
## asim_hirs4_n19_20181120190000.ensmean' has size 0. Returning a NULL data.table.
## [1] "analisis/diagfiles/E7/asim_hirs4_n19_20181120200000.ensmean"
## [1] "analisis/diagfiles/E7/asim_hirs4_n19_20181120210000.ensmean"
## Warning in fread(f): File 'analisis/diagfiles/E7/
## asim_hirs4_n19_20181120210000.ensmean' has size 0. Returning a NULL data.table.
## [1] "analisis/diagfiles/E7/asim_hirs4_n19_20181120220000.ensmean"
## [1] "analisis/diagfiles/E7/asim_hirs4_n19_20181120230000.ensmean"
## Warning in fread(f): File 'analisis/diagfiles/E7/
## asim_hirs4_n19_20181120230000.ensmean' has size 0. Returning a NULL data.table.
## [1] "analisis/diagfiles/E7/asim_hirs4_n19_20181121000000.ensmean"
## Warning in fread(f): File 'analisis/diagfiles/E7/
## asim_hirs4_n19_20181121000000.ensmean' has size 0. Returning a NULL data.table.
## [1] "analisis/diagfiles/E7/asim_hirs4_n19_20181121080000.ensmean"
## [1] "analisis/diagfiles/E7/asim_hirs4_n19_20181121200000.ensmean"
## [1] "analisis/diagfiles/E7/asim_hirs4_n19_20181121220000.ensmean"
# cat("Archivo ", basename(files[f]))

colnames(diag) <- c("sensor", "channel", "lat", "lon", "press", "dhr", "tb_obs", "tbc", "tbcnob",
                    "errinv", "qc", "emis", "tlapchn", "rzen", "razi", "rlnd", "rice", "rsnw", "rcld", 
                    "rcldp", paste0("pred", seq(8)), "date")

diag[sensor == "amsua_n15" & channel %in% c(1, 2, 3, 6, 15)] %>% 
  .[date == date[1]] %>% 
  ggplot(aes(ConvertLongitude(lon), lat)) +
  geom_point(aes(color = tbc - tbcnob)) +
  scale_color_divergent() +
  geom_mapa() +
  coord_sf(xlim = c(-75, -52), ylim = c(-42,-20)) +
  facet_wrap(~channel, ncol = 3) +
  theme_minimal()

diag[sensor == "amsua_n15" & channel == 2] %>% 
  .[date == date[1]] %>% 
  melt(measure.vars = paste0("pred", seq_len(8))) %>% 
  ggplot(aes(ConvertLongitude(lon), lat)) +
  geom_point(aes(color = value)) +
  scale_color_viridis_c() +
  geom_mapa() +
  coord_sf(xlim = c(-75, -52), ylim = c(-42,-20)) +
  facet_wrap(~variable, ncol = 4) +
  theme_minimal()

diag[sensor %in% c("amsua_aqua", "amsua_metop-a", "amsua_n15", "amsua_n18")] %>% 
  ggplot(aes(factor(channel), tbc - tbcnob)) +
  geom_boxplot() +
  facet_wrap(~sensor)

# diag[sensor %in% c("airs_aqua") & channel < 300] %>% 
#   ggplot(aes(channel, tbc - tbcnob)) +
#   geom_boxplot(aes(group = channel))